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ABSTRACT 

We propose a model to study the plasticity of ice single crystals by numerical simulations. 
The model includes the long-range character of the interaction among dislocations, as well 
as the possibility of mutual annihilation of these line defects characterized by its Burgers 
vector. A multiplication mechanism representing the activation of Frank- Read sources due 
to dislocation pinning is also introduced in the model. 

With our approach we are able to probe the dislocation patterns, which result from the 
dislocation dynamics. Furthermore, our results exhibit features characteristic of driven dy- 
namic critical phenomena such as scaling behavior, and avalanche dynamics. Some of these 
results account for the experimental findings reported for ice single crystals under creep 
deformation, like the power-law distributions of the acoustic emission intensity observed 
sistematically in experiments. 

INTRODUCTION 

The viscoplastic deformation of crystalline materials, such as ice single crystals, involves 
the motion of a large number of dislocations. Although the dynamics of an individual 
dislocation is a fairly well characterized phenomenon [p], |2|, the collective behavior of a 
large number of these defects appears to be an amazingly rich but poorly understood 
problem. The interaction between a pair of dislocations can be attractive or repulsive, 
depending on the orientation of their respective Burgers vectors; it grows logarithmically 
with the interline distance and allows mutual annihilation of defects. Dislocations may be 
incorporated into a crystal in the growth process, affecting the topology of the whole lattice. 
Moreover, under deformation conditions, dislocations can penetrate into the material from 
the sample surfaces, or be generated by various mechanism, as for example, in what is 
usually called a Frank-Read source, which is activated after the pinning of a dislocation 
loop. 

When a material is deformed under constant load (creep experiment), and dislocation 
motion is the dominant mechanism for viscoplastic deformation (other possible source is for 
instance crack nucleation and propagation), a constant strain-rate regime usually follows 
after the initial transient stage. Orowan's relation 7 = p m bv is known to prevail under such 
conditions, where 7 is the strain of the sample, p m is the density of mobile dislocations, 



b is the Burgers' vector, and v is the mean velocity of the dislocations. Obviously, this is 
a mean-field relation which neglects temporal and spatial fluctuations of both the density 
and the velocity fields. As a result of their interactions, however, dislocations tend to move 
cooperatively giving rise to a rather complex and heterogeneous slip process. Dislocations 
move in groups to form slip bands. Moving dislocations can pile-up against stable dislo- 
cation configurations such as walls or boundaries, which may eventually break apart. In 
this process, fluctuations in the dislocation density and velocity may be comparable to or 
greater than the mean values, and consequently, of great importance. 

The complex character of the collective dislocation dynamics reveals itself in experi- 
ments of acoustic emission (AE) ||. Sudden local changes of inelastic strain generate AE 
waves. Weiss and Grasso [||] soon realised that ice single crystals are particularly well suited 
for the study of dislocation dynamics. The perfect transparency of this material, easily (by 
eye) allows to rule out the possibility of cracks being present in the material which will, 
otherwise, interfere with the dislocations motion. Given the amplitude threshold and the 
frequency range accesible to the experimental apparatus, the AE signals detected seem to 
correspond to the synchronous motion of several dislocations, likely to occur for example 
during the breakaway of a pile of these defects, or the activation of a multiplication source. 
The AE experiments, however, have only access to information resulting from the inter- 
play of various magnitudes. Thus, the physical interpretation of the generated AE waves 
remains a major difficulty, and consequently, constitutes the main motivation of our work. 

The proportionality between global AE activity and global strain-rate (Orowan's rela- 
tion) has been tested experimentally. More locally, the AE intensity is thought to depend 
on the number of dislocations involved in a plastic instability, their length, and velocity. 
Various measurements of the acoustic activity recorded during a stress-constant step show 
that the AE signal takes place in the form of bursts spanning a wide range of amplitude 
values. In particular, the AE amplitude is distributed according to a power law. This be- 
havior is a consequence of the collective motion of dislocations which spontaneously gives 
rise to an avalanche- like dynamics, typical of slowly driven dissipative systems. The power 
law distribution provides evidence of scale-free cooperative behavior, whose origin could be 
ascribed to nonequilibrium continuous phase transitions 0] or self-organized criticality ||. 

DESCRIPTION OF THE MODEL 

To characterize the plastic deformation of a material from the perspective of nonequilibrium 
statistical mechanics, we propose a simplified dynamic model to simulate the system in 
rather general conditions. 

Ice single crystals deform essentially by slip on the basal plane (0001) (that we call the 
xy plane), i.e. the motion of dislocation lines or loops takes place by gliding on the xy 
planes. The possible Burgers vectors b are the three lattice vectors of an hexagonal lattice. 
For the sake of simplicity, we study a two-dimensional model representing a cross section 
of the crystal which is perpendicular to the basal planes and parallel to one of the lattice 
vectors, that is, for example, the xz plane. In this way, the dislocations constrained to 
move in this plane have all Burgers vectors parallel to the chosen lattice vector b = (b, 0, 0) 
and move along fixed lines parallel to the x axis. We also consider that all N dislocations 
are of edge type, and that, on average over many realizations, the number of dislocations 
with positive and negative Burgers vectors is the same. 

Several simplified models containing similar basic ingredients have been proposed in the 



literature in the last few years || |7|, |8|, [ITJ. A basic feature common to most models, is 
that dislocations interact with each other through the long-range elastic stress field they 
produce in the host material. An edge dislocation with Burgers vector b located at the 
origin gives rise to a shear stress a s at a point (x, z) of the form 



x 2 + z 2 ) 2 ' 



where D = u/27r(l — a) is a coefficient involving the shear modulus u and the Poisson's 
ration a for the material. In our model, we further assume that the dislocation velocities are 
linearly proportional to the local stress. Experimental evidence supports such a relationship 
for low stress conditions, which is indeed the case in our model. Accordingly, the velocity 
of the nth dislocation, if an external shear stress a e is also applied, is given by 

Vn = b n (Y,< m -° e )- (2) 

As the number of dislocations in any real crystals exceeds by far the number of defects 
we can handle in a computer, one usually introduces periodic boundary conditions (PBC) 
to effectively extend the size of our system. Due to the long range character of the force ([!]), 
we have exactly evaluated the Ewald sums of this expression to account for the interaction 
of a dislocation with all the infinite periodic replicas of all the other dislocations in a finite 
box of dimensions L x L. Contrary to what is stated in some references 0, we do not find 
any spureous results coming from the implementation of PBC in this fashion. Instead, we 
obtain artificial results when using the "nearest image" approximation and the truncation 
that this approximation implies. 

When the distance between two dislocations is of the order of a few Burgers vectors, 
the high stress and strain conditions close to the dislocation core invalidate the results 
obtained from a linear elasticity theory (i.e. Eq .(jl])). In these instances, phenomenological 
nonlinear reactions describe more accurately the real behavior of dislocations in a crystal. 
In particular in our model, we account for the annihilation of dislocations with opposite 
Burgers vectors when the distance between them is shorter than 2b. Thus the core of one 
dislocation in our model has a radius of size b. 

Another important feature of any computer model is the implementation of a mechanism 
for the multiplication of dislocations in the sample. It is widely believed that the Frank- 
Read mechanism |IJ is the most relevant for a gliding process of dislocations under creep 
deformation. Indeed Frank-Read sources (FRS) have been observed in ice. In a FRS 
multiplication occurs by pinning of a dislocation segment on the basal planes due, for 
example, to a defect in the crystal, or to dislocation dipoles, piles, and walls. Under an 
applied stress, the pinned segment bows out by glide and, if the local stress concentration 
is less than a critical value, a metastable configuration is attained where the line tension 
balances the stress. Beyond this threshold value, the dislocation segment wraps around 
itself, creating a new dislocation loop and restoring the original configuration. Thus a 
sequence of loops forms continuously from the source until the local shear stress drops 
below the activation value. In our model, we simulate this mechanism phenomenologically: 
a dislocation pair is generated (i) when the fraction of immobile dislocations is high (pinning 
is then more likely) and (ii) the local stress is large compared to a threshold value. Rather 
than fixing particular values of the parameters, we use a probabilistic procedure, keeping in 
mind that the details of the rule should not change the collective properties of the system. 
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Figure 1: a) A stable arrangement of dislocations in a box of size L = 3006. Various structures 
like dislocation dipoles, piles, and dislocation walls are present, b) The corresponding elastic 
stress map (absolute value) in a graylevel scale. Dark color indicates a low-stress region and 
lighter colors areas of higher stress. 

Annihilation and multiplication processes imply that the number of dislocations N in 
our system is not fixed in the course of time. Starting from a random configuration of 
dislocations, we let them relax until they find a stable arrangement (which could be an 
equilibrium configuration, or a long-lived metastable state). We solve numerically the A" 
equations of motion using an adaptive step size fifth-order Runge-Kutta method. So far, we 
have considered three different box sizes L = 1006, L = 2006, and L = 3006, with an initial 
number of dislocations A" = 400, A" = 800, and A" = 1500, respectively. After the system 
has reached a stable arrangement, the volume fraction of dislocations <fi = Nnb 2 /(Lb) 2 
ranges between 1 — 10%. Once in these conditions, we apply a small external shear stress 
and keep track of the various quantities describing the dynamics of the dislocations. 

RESULTS AND DISCUSSION 

Persistent slip bands or cell structures appearing on various length scales, are just a few 
examples of dislocation patterns. With our approach we are able to recover some of these 
structures. This is a subject which has raised a considerable amount of activity recently 
(see for example Refs. |0, [12] and references therein). Fig. |I]a) represents a common 
configuration of dislocations in stable conditions. There, we can observe the formation 
of several dipoles, piles, as well as dislocation walls delimiting various slip bands. The 
corresponding elastic stress map is depicted in Fig. [L]b). We plot the absolute value of 
the shear stress in a graylevel scale. The dark areas represent low-stress regions, and are 
consequently predominant in a stable configuration. The light portions mark a few regions 
of higher stress concentration given this particular arrangement. 

Second, we have been able to keep track of several quantities which play a key role in 
promoting scale invariant behavior, like the root-mean-square velocity, the local stress, the 
average number of dislocations, etc. In Figure [|a), we represent, for example, the root- 
mean-square velocity V m of a single run of the creep process. A curious feature is that after 
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Figure 2: a) Mean velocity as a function of time in our model of collective dislocation dynamics, 
b) The cumulative avalanche size distribution for three system sizes L = 100, 200, 300 represented 
in a double logarithmic scale. 



a burst, the relaxation of the velocity V m to a still configuration is power-law like, i.e. slow 
and without any characteristic time. 

We have defined the size s of an acoustic avalanche as the sum of the root-mean-square 
velocity of all the moving dislocations, that is s = Yl \ v i\- The cumulative distribution P c (s) 
of the avalanches obtained after averaging over several realizations is depicted in Figure 0b) 
for the three system sizes studied L = 100, 200, 300. We recover a very clear power law 
distribution (P c (s) ~ s _r ) extending over close to two decades. The exponent r ~ 0.65 
is in reasonable agreement with experimental data ||. The distribution cut-off for large 
values of s is due to the finite size of the sample. (As one would expect, the cut-off is 
scaling accordingly to the sytem size.) This clearly points out the presence of a very large 
(or infinite) characteristic size for the acoustic events. It is worth remarking that large 
stresses introduce a characteristic scale in the process. A more detailed study as a function 
of the applied stress is in progress [[U| . 

The numerical investigation of the present model provides striking evidences for the 
collective critical behavior of dislocation motion under external stress. Relevant magni- 
tudes characterizing the dynamics of dislocations show power law behavior signalling the 
absence of any characteristic length or time in the process. The response to an infinitesimal 
perturbation (slow injection of new dislocations) exhibits singular behavior in the guise of 
avalanches distributed over many length scales. Avalanche dynamics is the rule rather 
the exception in slowly driven disordered systems. Examples can be found in fracturing 
of wood and concrete ||13||, Barkhausen effect and flux lines in high-T c superconduc- 



tors [ i~5"| . Under the external drive (the stress in the present case), the system jumps 
between metastable or pinned configurations in which the dynamics is virtually frozen. In 
the limit of a very slow driving (that is equivalent to a fine tuning close to the depinning 
critical point), the disordered energy landscape is explored quasistatically and the response 
function exhibits critical properties ||16|| . Tipically, a basic ingredient for this behavior is 
the presence of quenched disorder acting as the source of pinning in the system. Noticeably, 
the system under study does not contain any external source of disorder. Pinned states are 
due to the various dislocation patterns that spontaneously develop in the system. Struc- 



tures such as dipoles, piles, and dislocation walls, play the role of self-generated pinning 
centers that create the pinning force landscape. The new scenario poses many new and 
interesting questions for a definitive identification and understanding of the critical nature 
of dislocation dynamics. 

ACKNOWLEDGMENTS 

We gratefully acknowledge R. Pastor-Satorras advice on the algorithm implementation. We 
also thank J.R. Grasso and J.Weiss for fruitful discussions regarding their experiment, and 
for providing us with the experimental data prior to publication. 

References 

[1] J. P. Hirth and J. Lothe, Theory of Dislocations (Krieger Publishing Company, 1992). 
[2] F.R.N. Nabarro, Theory of Crystal Dislocations (Dover, New York, 1992). 
[3] J. Weiss and J.R. Grasso, J. Phys. Chem. B 101, 6113 (1997). 

[4] Phase Transitions and Critical Phenomena, edited by C. Domb and M.S. Green (Aca- 
demic Press, London 1972-1976), Vols. 1-17. 

[5] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987); A. Vespignani 
and S.Zapperi, Phys. Rev. E, 57, 6345 (1998). 

[6] J. Lepinoux and L.P. Kubin, Scripta Metall. 21, 833 (1987). 

[7] R.J. Amodeo and N.M. Ghoniem, Phys. Rev. B 41, 6958 and 6968 (1990). 

[8] I. Groma and G.S. Pawley, Phil. Mag. A, 67, 1459 (1993). 

[9] R. Fournet and J.M. Salazar, Phys. Rev. B 53, 6283 (1996). 
[10] M.C. Miguel, A. Vespignani, and S. Zapperi, in preparation. 
[11] P. Hahner, K. Bay, and M. Zaiser, Phys. Rev. Lett. 81, 2470 (1998). 
[12] B. Bako and I. Groma, Phys. Rev. B 60, 122 (1999). 

[13] A. Garcimartin, A. Guarino, L. Bellon, and S. Ciliberto, Phys. Rev. Lett. 79, 3202 
(1997); A. Petri, G. Paparo, A. Vespignani, A. Alippi, and M. Costantini Phys. Rev. 
Lett. 73, 3423 (1994). 

[14] G. Bertotti, G. Durin, and A. Magni, J. Appl. Phys. 75, 5490 (1994). 

[15] S. Field, J. Witt, F. Nori and X. Ling. Phys. Rev. Lett. 74, 1206 (1995). 

[16] For a review on mechanisms that generate spontaneously critical behavior see R. Dick- 
man, M. A. Munoz, A. Vespignani and S. Zapperi, e-print |cond-mat/ 99 10454 . 



